#include <math.h>

#include "fdecl.h"

void compute_element_stiffness_matrix(hybroid *h)
{
  int i, j;
  double ww1 = h->w1, ww0 = h->w0, ww04 = 4.0 * ww0;
  double ust = h->ustep, vst = h->vstep;
  double ust2 = ust * ust, vst2 = vst * vst;
  double area3 = ust * vst / 2.0, area4 = ust * vst / 4.0; 


  /* compute upper diagonal portion of the 4 node stiffness matrix */

  h->s400 = area4 * 4.0 * (ww1 / (3.0 * ust2) +
			       ww1 / (3.0 * vst2) + ww0);
  h->s401 = area4 * (-4.0 * ww1 / (3.0 * ust2) +
			 2.0 * ww1 / (3.0 * vst2) + ww04);
  
  h->s402 = area4 * (- 2.0 * ww1 / (3.0 * ust2) -
			 2.0 * ww1 / (3.0 * vst2) + ww04);
  h->s403 = area4 * (2.0 * ww1 / ( 3.0 * ust2) -
			 4.0 * ww1 / ( 3.0 * vst2) + ww04);
  h->s411 = area4 * 4.0 * (ww1 / (3.0 * ust2) +
			       ww1 / (3.0 * vst2) + ww0);
  h->s412 = area4 * (2.0 * ww1 / (3.0 * ust2) -
			 4.0 * ww1 / (3.0 * vst2) + ww04);
  h->s413 = area4 * (-2.0 * ww1 / (3.0 * ust2) - 
			 2.0 * ww1 / (3.0 * vst2) + ww04);
  h->s422 = area4 * 4.0 * (ww1 / (3.0 * ust2) + 
			       ww1 / (3.0 * vst2) + ww0);
  h->s423 = area4 * (-4.0 * ww1 / (3.0 * ust2) +
			 2.0 * ww1 / (3.0 * vst2) + ww04);
  h->s433 = area4 * 4.0 * (ww1 / (3.0 * ust2) + 
			       ww1 / (3.0 * vst2) + ww0);
  
  /* compute upper diagonal portion of the 3 node stiffness matrix */

  // south pole elements
  h->s300_1 = area3 * (ww1 / ust2 + ww0);
  h->s301_1 = area3 * (-ww1 / ust2 + ww0);
  h->s302_1 = area3 * ww0;
  h->s311_1 = area3 * (ww1 / ust2 + ww1 / vst2 + ww0);
  h->s312_1 = area3 * (-ww1 / vst2 + ww0);
  h->s322_1 = area3 * (ww1 / vst2 + ww0);

  // north pole elements
  h->s300_2 = area3 * (ww1 / ust2 + ww1 / vst2 + ww0);
  h->s301_2 = area3 * (-ww1 / ust2 + ww0);
  h->s302_2 = area3 * (-ww1 / vst2 + ww0);
  h->s311_2 = area3 * (ww1 / ust2 + ww0);
  h->s312_2 = area3 * ww0;
  h->s322_2 = area3 * (ww1 / vst2 + ww0);

  // elements in the mesh
  h->s300_3 = area3 * (ww1 / vst2 + ww0);
  h->s301_3 = area3 * (-ww1 / vst2 + ww0);
  h->s302_3 = area3 * ww0;
  h->s311_3 = area3 * (ww1 / ust2 + ww1 / vst2 + ww0);
  h->s312_3 = area3 * (-ww1 / ust2 + ww0);
  h->s322_3 = area3 * (ww1 / ust2 + ww0);
  
}

void compute_element_shape_functions_at_integration_points(hybroid *h)
{

/* compute Gauss point values for the shape functions of 4 node element */

  h->n4_1[0] = 0.25 * (1.0 + 1.0 / sqrt(3.0)) * (1.0 + 1.0 / sqrt(3.0));
  h->n4_1[1] = 0.25 * (1.0 + 1.0 / sqrt(3.0)) * (1.0 - 1.0 / sqrt(3.0));
  h->n4_1[2] = 0.25 * (1.0 - 1.0 / sqrt(3.0)) * (1.0 + 1.0 / sqrt(3.0));
  h->n4_1[3] = 0.25 * (1.0 - 1.0 / sqrt(3.0)) * (1.0 - 1.0 / sqrt(3.0));

  h->n4_2[0] = 0.25 * (1.0 - 1.0 / sqrt(3.0)) * (1.0 + 1.0 / sqrt(3.0));
  h->n4_2[1] = 0.25 * (1.0 - 1.0 / sqrt(3.0)) * (1.0 - 1.0 / sqrt(3.0));
  h->n4_2[2] = 0.25 * (1.0 + 1.0 / sqrt(3.0)) * (1.0 + 1.0 / sqrt(3.0));
  h->n4_2[3] = 0.25 * (1.0 + 1.0 / sqrt(3.0)) * (1.0 - 1.0 / sqrt(3.0));

  h->n4_3[0] = 0.25 * (1.0 - 1.0 / sqrt(3.0)) * (1.0 - 1.0 / sqrt(3.0));
  h->n4_3[1] = 0.25 * (1.0 - 1.0 / sqrt(3.0)) * (1.0 + 1.0 / sqrt(3.0));
  h->n4_3[2] = 0.25 * (1.0 + 1.0 / sqrt(3.0)) * (1.0 - 1.0 / sqrt(3.0));
  h->n4_3[3] = 0.25 * (1.0 + 1.0 / sqrt(3.0)) * (1.0 + 1.0 / sqrt(3.0));

  h->n4_4[0] = 0.25 * (1.0 + 1.0 / sqrt(3.0)) * (1.0 - 1.0 / sqrt(3.0));
  h->n4_4[1] = 0.25 * (1.0 + 1.0 / sqrt(3.0)) * (1.0 + 1.0 / sqrt(3.0));
  h->n4_4[2] = 0.25 * (1.0 - 1.0 / sqrt(3.0)) * (1.0 - 1.0 / sqrt(3.0));
  h->n4_4[3] = 0.25 * (1.0 - 1.0 / sqrt(3.0)) * (1.0 + 1.0 / sqrt(3.0));

/* compute Radau point values for the shape functions of 3 node element */

/* south pole elements */

  h->n31_1[0] = 0.5;
  h->n31_1[1] = 0.0;
  h->n31_1[2] = 0.5;

  h->n31_2[0] = 0.5;
  h->n31_2[1] = 0.5;
  h->n31_2[2] = 0.0;

  h->n31_3[0] = 0.0;
  h->n31_3[1] = 0.5;
  h->n31_3[2] = 0.5;
  
/* north pole elements */

  h->n32_1[0] = 0.5;
  h->n32_1[1] = 0.0;
  h->n32_1[2] = 0.5;

  h->n32_2[0] = 0.5;
  h->n32_2[1] = 0.5;
  h->n32_2[2] = 0.0;

  h->n32_3[0] = 0.0;
  h->n32_3[1] = 0.5;
  h->n32_3[2] = 0.5;
}










